Average wage in Russia


In [1]:
from __future__ import division

import numpy as np
import pandas as pd

from scipy import stats
import statsmodels.api as sm

%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns

import warnings
from itertools import product
from datetime import *
from dateutil.relativedelta import *

from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"


C:\ProgramData\Anaconda3\envs\python2\lib\site-packages\statsmodels\compat\pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.
  from pandas.core import datetools

Для выполнения этого задания нам понадобятся данные о среднемесячных уровнях заработной платы в России:


In [2]:
#Reading data
wage = pd.read_csv('WAG_C_M.csv', sep=';', index_col='month', parse_dates=True, dayfirst=True)
wage.info()


<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 294 entries, 1993-01-01 to 2017-06-01
Data columns (total 1 columns):
WAG_C_M    294 non-null float64
dtypes: float64(1)
memory usage: 4.6 KB

In [3]:
wage.head()


Out[3]:
WAG_C_M
month
1993-01-01 15.3
1993-02-01 19.1
1993-03-01 23.6
1993-04-01 30.6
1993-05-01 37.5

In [4]:
_ = plt.figure(figsize=(15,7))
_ = wage.WAG_C_M.plot()
_ = plt.title('Average nominal wage')


Проверка стационарности и STL-декомпозиция ряда:


In [5]:
_ = sm.tsa.seasonal_decompose(wage.WAG_C_M).plot()
print('Augmented Dickey-Fuller unit root test p=%f' % sm.tsa.stattools.adfuller(wage.WAG_C_M)[1])


Augmented Dickey-Fuller unit root test p=0.995929

Стабилизация дисперсии

Сделаем преобразование Бокса-Кокса для стабилизации дисперсии:


In [6]:
wage['WAG_C_M_box'], lmbda = stats.boxcox(wage.WAG_C_M)

_ = plt.figure(figsize=(15,7))
_ = wage.WAG_C_M_box.plot()
_ = plt.title(u'Transformed average nominal wage')

print('Optimal parameter of the Box-Cox power transformation: %f' % lmbda)
print('Augmented Dickey-Fuller unit root test p=%f' % sm.tsa.stattools.adfuller(wage.WAG_C_M_box)[1])


Optimal parameter of the Box-Cox power transformation: 0.270218
Augmented Dickey-Fuller unit root test p=0.719046

Стационарность

Критерий Дики-Фуллера отвергает гипотезу нестационарности, но визуально в данных виден тренд. Попробуем сезонное дифференцирование; сделаем на продифференцированном ряде STL-декомпозицию и проверим стационарность:


In [7]:
wage['WAG_C_M_box_diff'] = wage.WAG_C_M_box - wage.WAG_C_M_box.shift(12)
_ = sm.tsa.seasonal_decompose(wage.WAG_C_M_box_diff.dropna()).plot()
print('Augmented Dickey-Fuller unit root test p=%f' % sm.tsa.stattools.adfuller(wage.WAG_C_M_box_diff.dropna())[1])


Augmented Dickey-Fuller unit root test p=0.009562

Критерий Дики-Фуллера отвергает гипотезу нестационарности, НО полностью избавиться от тренда не удалось. Попробуем добавить ещё обычное дифференцирование:


In [8]:
wage['WAG_C_M_box_diff2'] = wage.WAG_C_M_box_diff - wage.WAG_C_M_box_diff.shift(1)
_ = sm.tsa.seasonal_decompose(wage.WAG_C_M_box_diff2.dropna()).plot()
print('Augmented Dickey-Fuller unit root test p=%f' % sm.tsa.stattools.adfuller(wage.WAG_C_M_box_diff2.dropna())[1])


Augmented Dickey-Fuller unit root test p=0.000000

Гипотеза нестационарности отвергается с ещё большим уровнем значимости, и визуально ряд выглядит лучше — тренда больше нет.

Подбор модели

Посмотрим на ACF и PACF полученного ряда:


In [9]:
plt.figure(figsize=(15,10))

ax = plt.subplot(211)
sm.graphics.tsa.plot_acf(wage.WAG_C_M_box_diff2.dropna()[12:].squeeze(), lags=50, ax=ax);

ax = plt.subplot(212)
sm.graphics.tsa.plot_pacf(wage.WAG_C_M_box_diff2.dropna()[12:].squeeze(), lags=50, ax=ax);


Начальные приближения: Q=0, q=1, P=1, p=1.


In [10]:
ps = range(0, 2)
d=1
qs = range(0, 2)
Ps = range(0, 2)
D=1
Qs = range(0, 1)

In [11]:
parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
parameters_list
len(parameters_list)


Out[11]:
[(0, 0, 0, 0),
 (0, 0, 1, 0),
 (0, 1, 0, 0),
 (0, 1, 1, 0),
 (1, 0, 0, 0),
 (1, 0, 1, 0),
 (1, 1, 0, 0),
 (1, 1, 1, 0)]
Out[11]:
8

In [12]:
%%time
results = []
best_aic = float("inf")
warnings.filterwarnings('ignore')

for param in parameters_list:
    #try except нужен, потому что на некоторых наборах параметров модель не обучается
    try:
        model=sm.tsa.statespace.SARIMAX(wage.WAG_C_M_box, order=(param[0], d, param[1]), 
                                        seasonal_order=(param[2], D, param[3], 12)).fit(disp=-1)
    #выводим параметры, на которых модель не обучается и переходим к следующему набору
    except ValueError:
        print('wrong parameters:', param)
        continue
    aic = model.aic
    #сохраняем лучшую модель, aic, параметры
    if aic < best_aic:
        best_model = model
        best_aic = aic
        best_param = param
    results.append([param, model.aic])
    
warnings.filterwarnings('default')


('wrong parameters:', (0, 0, 0, 0))
Wall time: 1.21 s

In [13]:
result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
print(result_table.sort_values(by = 'aic', ascending=True).head())


     parameters        aic
4  (1, 0, 1, 0)  36.753491
2  (0, 1, 1, 0)  37.356642
6  (1, 1, 1, 0)  38.745133
3  (1, 0, 0, 0)  40.310501
1  (0, 1, 0, 0)  40.901350

Лучшая модель:


In [14]:
print(best_model.summary())


                                 Statespace Model Results                                 
==========================================================================================
Dep. Variable:                        WAG_C_M_box   No. Observations:                  294
Model:             SARIMAX(1, 1, 0)x(1, 1, 0, 12)   Log Likelihood                 -15.377
Date:                            Thu, 17 Aug 2017   AIC                             36.753
Time:                                    23:55:41   BIC                             47.804
Sample:                                01-01-1993   HQIC                            41.179
                                     - 06-01-2017                                         
Covariance Type:                              opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1         -0.1643      0.045     -3.633      0.000      -0.253      -0.076
ar.S.L12      -0.1471      0.045     -3.238      0.001      -0.236      -0.058
sigma2         0.0653      0.004     15.288      0.000       0.057       0.074
===================================================================================
Ljung-Box (Q):                       50.43   Jarque-Bera (JB):                50.74
Prob(Q):                              0.12   Prob(JB):                         0.00
Heteroskedasticity (H):               1.55   Skew:                             0.14
Prob(H) (two-sided):                  0.04   Kurtosis:                         5.06
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

Её остатки:


In [15]:
_ = plt.figure(figsize=(15,12))
_ = plt.subplot(211)
_ = best_model.resid[13:].plot()
_ = plt.ylabel(u'Residuals')

_ = ax = plt.subplot(212)
_ = sm.graphics.tsa.plot_acf(best_model.resid.values.squeeze(), lags=50, ax=ax)

print("Критерий Стьюдента: p=%f" % stats.ttest_1samp(best_model.resid[13:], 0)[1])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(best_model.resid[13:])[1])


Критерий Стьюдента: p=0.136758
Критерий Дики-Фуллера: p=0.000070

Остатки несмещены (подтверждается критерием Стьюдента), стационарны (подтверждается критерием Дики-Фуллера и визуально), неавтокоррелированы (подтверждается критерием Льюнга-Бокса и коррелограммой). Посмотрим, насколько хорошо модель описывает данные:


In [16]:
def invboxcox(y,lmbda):
    if lmbda == 0:
        return(np.exp(y))
    else:
        return(np.exp(np.log(lmbda*y+1)/lmbda))

In [17]:
wage['model'] = invboxcox(best_model.fittedvalues, lmbda)

_ = plt.figure(figsize=(15,7))
_ = wage.WAG_C_M.plot()
_ = wage.model[13:].plot(color='r')
_ = plt.title('Average nominal wage')


Прогноз


In [18]:
wage2 = wage[['WAG_C_M']]
date_list = [datetime.strptime("2017-07-01", "%Y-%m-%d") + relativedelta(months=x) for x in range(0,36)]
future = pd.DataFrame(index=date_list, columns=wage2.columns)
wage2 = pd.concat([wage2, future])
wage2['forecast'] = invboxcox(best_model.predict(start=294, end=329), lmbda)

_ = plt.figure(figsize=(15,7))
_ = wage2.WAG_C_M.plot()
_ = wage2.forecast.plot(color='r')
_ = plt.title('Average nominal wage')